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We derive the equations of motion, the periastron shift, and the gravitational radiation damping 
for quasicircular compact binaries in a massive variant of the Brans-Dicke theory of gravity. We also 
study the Shapiro time delay and the Nordtvedt effect in this theory. By comparing with recent 
observational data, we put bounds on the two parameters of the theory: the Brans-Dicke coupling 
parameter lobu and the scalar mass m^- We find that the most stringent bounds come from Cassini 
measurements of the Shapiro time delay in the Solar System, that yield a lower bound wbd > 40000 
for scalar masses rus < 2.5 x 10~^''eV (or Compton wavelengths As = h/{msc) > 5 x 10^° km), to 
95% confidence. In comparison, observations of the Nordtvedt effect using Lunar Laser Ranging 
(LLR) experiments yield cjbd > 1000 for rUs < 2.5 x 10~^''eV. Observations of the orbital period 
derivative of the quasicircular white dwarf-neutron star binary PSR J1012-I-5307 yield uj-bo > 1250 
for rUs < lO^^^eV (As > 1.2 x 10^^ km). A first estimate suggests that bounds comparable to 
the Shapiro time delay may come from observations of radiation damping in the eccentric white 
dwarf-neutron star binary PSR J1141-6545, but a quantitative prediction requires the extension of 
our work to eccentric orbits. 



General relativity (GR) occupies a well earned place 
next to the standard model as one of the two pillars of 
modern physics. All observational evidence to date sup- 
ports GR as the correct classical theory of gravitation, 
but there are countless attempts at developing alterna- 
tive theories of gravity. Two of the main motivations 
for these efforts are the desire to formulate a fully quan- 
tizable theory of gravity, and the quest to uncover the 
mechanisms underlying the dark energy problem in cos- 
mology. In addition, the vast majority of tests of GR that 
have been carried out to date are in the weak-field, low 
energy regime, but it is widely believed that GR may in- 
deed break down at higher energies. The direct observa- 
tion of gravitational waves with Earth- and space-based 
detectors will mark the dawn of a new era, allowing us 
to probe gravity in the dynamical, strong-field regime. 
For these reasons, the study of gravitational radiation in 
modified theories of gravity has become a central issue. 

One of the most popular and simple alternative the- 
ories of gravity is scalar-tensor theory, in which gravity 
is mediated by both a scalar and a tensor field, coupled 
together in a nontrivial manner through the presence of 
a nonminimal coupling term in the action [lH3|. The 
existence of scalar partners to the graviton is predicted 
in all extra-dimensional theories, and scalar fields play a 
crucial role in modern cosmology. Scalar-tensor theories 
are consistent, have a well-posed Cauchy problem, and 
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respect many of the symmetries of GR. They are also 
conformally equivalent to GR (if the coupling with mat- 
ter is nonstandard), allowing us to employ the same tech- 
niques used to solve the Einstein field equations as long 
as we work in the Einstein frame Finally, generic 

scalar-tensor theories can be shown to be equivalent to 
f{R) theories i,!- A good account of the motivations 
behind scalar-tensor theories, including their historical 
development, can be found in [llQ. 

String theory suggests the existence of massive but 
light scalar fields ("axions") with masses possibly as 
small as the Hubble scale (~ 10"'^^ eV). If we do indeed 
live in a "string axiverse" , CMB observations, galaxy sur- 
veys and measurements of black hole spins may offer ex- 
citing experimental opportunities to set constraints on 
the mass of these scalars 0, 0] • 

Here we are interested in the possibility of constrain- 
ing the mass and coupling of massive scalars via present 
(electromagnetic) and future (gravitational-wave) obser- 
vations of compact binaries. Until recently, calculations 
of gravitational radiation damping in scalar-tensor theo- 
ries (see e.g. [8l-[T]|) have focused mostly on the massless 
case. Due to the interest of light scalars in cosmology and 
high-energy physics, this restriction has been dropped in 
more recent work. For example it has been shown that 
resonant, superradiant effects induced by light, massive 
scalars may produce "floating orbits" when small com- 
pact objects inspiral into rotating black holes, leaving a 
distinct signature in gravitational waves [l^, [13 ■ 

A commonly held belief is that only mixed binaries 
(i.e., binaries whose members have different gravitational 
binding energy) can produce significant amounts of scalar 
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gravitational radiation. There are two reasons for this. 
The first is that, under standard assumptions, dipole ra- 
diation is produced due to violations of the strong equiv- 
alence principle when the binary members have unequal 
"sensitivities": si ^ S2- These sensitivities are defined in 
Eq. ([Tl]) below, and they are related to the gravitational 
binding energy of each binary member. In other words, 
dipole radiation is produced when the system's center of 
mass is offset with respect to the center of inertia (see 
G-g- 3), so that mixed binaries and eccentric binaries 
would be the best target to constrain scalar-tensor theo- 
ries. The second reason is the black hole no-hair theorem, 
i.e. the fact that black hole solutions in scalar-tensor 
theories are the same as in GR (see [31 and references 
therein). Building on earlier work by Jacobson [l5l |. Hor- 
batsch and Burgess recently pointed out that slowly vary- 
ing scalar fields may violate the no-hair theorem, so that 
even black hole-black hole binaries may produce dipole 
radiation [l^. They also developed a formalism to test 
generic scalar-tensor theories using binary pulsars [13 ■ 

For all these reasons, a study of gravitational radiation 
in massive scalar-tensor theories is quite timely. In this 
paper we derive the period derivative due to scalar and 
tensor radiation in theories with a massive scalar field. 
For simplicity we focus on circular binaries, but (as we 
will see below) the generalization of our results to eccen- 
tric binaries would be of great observational interesl[3. 

For the reader's convenience, here we give an execu- 
tive summary of our main results. Consider a compact 
binary in circular orbit with component masses rrii and 
sensitivities Si (i = 1 , 2). Then the period derivative due 
to the emission of scalar and tensor gravitational waves 
in the massive Brans-Dicke theory is 
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Q is the Heaviside function, r is the separation of the 
binary members, rus is the mass of the scalar field, m = 
nil -\- 7712 and fi = mim2/m are the total and reduced 
masses of the system, S = S2 — si and furthermore 

^ 2 + cjbd ' 

G = 1 - i{si + S2 - 2siS2) , 

sim2 + miS2 
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^ We will be working in units h = c = G = 1 throughout the 
paper. Greek indices will span both spatial and time components 
0, 1, 2, 3. Roman indices run over the spatial components 1, 2, 3 
only. We will adopt the metric signature (— , -|-, -|-, -|-). 



Note that scalar dipole radiation is emitted only when 
the binary's orbital frequency uj > uis and the difference 
in sensitivities S ^ 0, while scalar quadrupole/monopole 
radiation is emitted only when 2uj > nis and it also van- 
ishes for two black holes (since in that case si = S2 = 1/2 
and r = 0). This result is only strictly valid in the limit of 
a very massive (m^r ^ 1) or very light (mgr <^ 1) scalar. 
However corrections due to an intermediate mass scalar 
always enter with at least a factor of the small parameter 
^, so this should be a relatively good approximation for 
the full range of scalar masses. 

In addition to deriving the orbital period derivative 
due to gravitational radiation, we also revisit the calcu- 
lations of the Shapiro time delay and of the Nordtvcdt 
effect in the massive Brans-Dicke theory. As we will 
see, the presence of the massive scalar does not allow 
a straightforward implementation of the parametrized 
post-Newtonian formalism. By comparing our results 
for the orbital period derivative, Shapiro time delay and 
Nordtvedt parameter against recent observational data, 
we put constraints on the parameters of the theory: the 
scalar mass TOs and the Brans-Dicke coupling parameter 
cjbd • Our bounds are summarized in Figure [TJ 

Wc find that the most stringent bounds come from the 
observations of the Shapiro time delay in the Solar Sys- 
tem provided by the Cassini mission (these bounds were 
already studied by Perivolaropoulos, although he used a 
slightly different notation [13 )■ From the Cassini obser- 
vations we obtain wbd > 40000 for to^ < 2.5 x 10~2°eV, 
to 95% confidence. Observations of the Nordtvedt ef- 
fect using the Lunar Laser Ranging (LLR) experiment 
yield a slightly weaker bound of wbd > 1000 for nis < 
2.5 X 10~2°eV. Observations of the orbital period deriva- 
tive of the circular white-dwarf neutron-star (WD-NS) 
binary system PSR J1G12+5307 yields wbd > 1250 for 
nis < 10~2°eV. The limiting factor here is our abil- 
ity to obtain precise measurements of the masses of the 
component stars as well as of the orbital period deriva- 
tive, once kinematic corrections have been accounted for. 
However, there is considerably more promise in the ec- 
centric binary system PSR J1141-6545. This system has 
allowed for remarkably precise measurements of the or- 
bital period derivative, of the component star masses and 
of the periastron shift, making it a promising candidate 
for constraining alternative theories of gravity. Unfortu- 
nately the system has nonnegligible eccentricity. Gen- 
eralizing our result for the orbital period derivative to 
eccentric binaries is a significant (but worthy) algebraic 
undertaking. 

The plan of the paper is as follows. In section U we 
describe and motivate the Brans-Dicke theory with a 
massive scalar field. In section |ll] we perform a post- 
Newtonian expansion of the field equations. In section 
mil we deal with the Shapiro time delay. In section IIVI 
we proceed to obtain the equations of motion of a binary 
system as well as the periastron shift. In section |V] we 
discuss the Nordtvedt effect. In scction lVTl we give details 
of the derivation of the gravitational radiation damping 
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FIG. 1. Left: Lower bound on (lobo + 3/2) as a function of the mass of the scalar rUs from the Cassini mission data (black 
solid line; cf. [3), period derivative observations of PSR J1141-6545 (dashed red line) and PSR J1012+5307 (dot-dashed green 
line), and Lunar Laser Ranging experiments (dotted blue line). Vertical lines indicate the masses corresponding to the typical 
radii of the systems: lAU (black solid line) and the orbital radii of the two binaries (dashed red and dot-dashed green lines). 
Right panel: upper bound on ^ as a function of ms. Linestyles are the same as in the left panel. Note that the theoretical 
bounds on the coupling parameters are lu > —3/2 and ^ < 2. 



of a compact binary system due to scalar and tensor grav- 
itational radiation. In section IVlIl we use these results to 
put bounds on the parameters of the theory. In the con- 
clusions we point out possible future extensions of our 
work. Appendix [A] outlines a step-by-step derivation of 
the post-Newtonian expansion of the scalar field and of 
the metric. Appendix [B] provides details on certain in- 
tegrals that appear in the calculation of the energy flux. 
Finally, Appendix [C] contains a short summary of com- 
pact binary observations relevant to this work. 



I. THE BRANS-DICKE THEORY WITH A 
MASSIVE SCALAR FIELD 

A. The generic scalar-tensor theory with a single 
scalar field 

A general class of scalar-tensor theories containing a 
single scalar field in addition to the tensor field was stud- 
ied by Bergmann and Wagoner [1, [l^ . We can charac- 
terize the Bergmann- Wagoner theory via the following 
postulates: 

1) The principle of general covariance is imposed, leading 
to tensorial equations. 

2) The field equations are derived from the action 



(3) 



where Jffc and are the Lagrangian densities for the 
gravitational and matter fields, respectively. 
3) We postulate that the long-range forces of nature are 
mediated by the three lowest spin bosons, and assume 
that the electromagnetic field is the only vector field. 



This leaves a scalar degree of freedom (f> and a tensor 
degree of freedom (the metric g^^) to describe the dy- 
namics of the gravitational field. 

4) The field equations are of at most of second differential 
order, and the tensor and scalar fields are nonminimally 
coupled; this leads us to the general form 

- i-g)^ [him + ;(0)5^>,;.0,. + m] (4) 

for the gravitational Lagrangian density, where h{(f>), l{4)) 
and A((/)) are arbitrary functions of the scalar field (f>. 

5) We postulate a principle of mutual coupling, in which 
the matter Lagrangian density depends on the gravita- 
tional fields according to 



(5) 



where 7/'((/)) is a fourth arbitrary function of 0, and ^ 
represents the collective matter fields. This guarantees 
consistency with the strong equivalence principle 

Now let us make the conformal transformation g^^, — >■ 
■0^(0)(7^i,, and in doing so move into a conformal frame 
in which the matter fields do not couple directly (but 
only indirectly, via the metric) to the scalar field: this is 
commonly referred to as the Jordan frame [l|, Q . Fur- 
thermore, without loss of generality we can redefine the 
scalar field such that h{(j)) — ^ cj). These two redefinitions 
recast the action into the form 



S 



1 

16^ 



(j>R- 



i-g)-^d'x 
(6) 



which has the additional advantage that the resulting 
weak-field equations for and decouple from one 



4 



another. The generic theory now contains two undeter- 
mined functions: the cosmological Junction Mid) and 
the coupling function ijo{(f)) (in the language of [2Ql). The 
effect of the couphng function on compact binary dy- 
namics has been studied extensively, and it can lead to 
interesting con seq uences if "spontaneous scalarization" 
occurs nil mH2J- Here we focus on the cosmological 
function, which has three major effects in the generic 
theory. Firstly, in the resulting field equations for g^j^ 
it plays the role of a cosmological constant. Secondly, it 
endows the scalar with mass: this manifests itself most 
clearly in the fact that solutions for (j) for isolated systems 
contain Yukawa-like terms e~"'='"/r, where is the mass 
of the scalar field, which in turn gives the field a char- 
acteristic range ^ ^ 1/ [l^l ■ Finally, the cosmological 
function may introduce nonlinearities in the dynamics of 
the scalar field. 



B. The matter action and the field equations 

Let us now turn to the matter action. Throughout this 
paper we will make the assumption that all bodies in our 
system can be treated as point masses. Einstein, Infeld 
and Hoffmann (EIH) [2^ developed a method for obtain- 
ing the equations of motion for a system of gravitating 
point-like masses. In their approach, one begins by ob- 
taining the local gravitational field of a single body (in a 
comoving frame) , under the assumption that the body is 
small and nearly spherical. One then proceeds to match 
the interbody gravitational fields onto the obtained lo- 
cal field of the single body under inspection; imposing 
self consistency yields the EIH equations of motion. The 
same equations of motion can be obtained with signif- 
icantly less effort, albeit at the sacrifice of some rigor, 
by taking the stress-energy tensor to be a distribution 
of delta functions and neglecting any infinite self-energy 
terms as they arise [l^l- In scalar-tensor theory, how- 
ever, we must deal with the additional complication that 
the inertial mass and internal structure of a gravitating 
body will depend on the local value of the scalar field 
(i.e. the local value of the effective gravitational "con- 
stant"). Variations in internal structure may act back 
on the motion of the body, leading to violations of the 
(weak) equivalence principle. Eardley [2^ showed that 
these effects could be accounted for by simply supposing 
that the masses of the bodies are in general functions of 

I 



,dM(£, 
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the scalar field, such that the matter action for a system 
of point-like masses can be written as 



5*1^1 = -^ J ma{(t))dTa, 



(7) 



where the particles (labeled by a) have inertial masses 
ma{(t>), and Tq is the proper time of particle a measured 
along its worldline x^. The distributional stress-energy 
tensor T''" and its trace T hence take the form 

T^'^ix') = i-g)-i J2 ^a{^)^S*{x^ - O , (8) 

a 

T = g.^T^'^ = -{-g)-h "^^^x' - <) • 0) 



Far from the system, the scalar will take on its cosmo- 
logically imposed value, denoted by ^o- The relationship 
between the effective gravitational constant, G, and the 
scalar field (j) is therefore (in our chosen system of units) 
G ~ (t'o/4'- III the post- Newtonian limit, we expand (j) 
about its asymptotic value and define the small pertur- 
bation (p such that = 00 + </3. In this case, we can write 
the variation of the inertial masses nia with as 



mai4>) =ma (In G) ma (00 
1 , 







(10) 



where we have defined the "first and second sensitivities" 
Sa and sjj to be0 



(9(ln TTla) 



d(ln G) 



92(ln TUa) 



(11) 



9(ln G)2 
The full action is now given by 

S~ I [</>i?-^5^''0,^0,. + M(0)](-5)^d^x 

-Y. [ T^a{<t>)dTa. (12) 

By varying the action with respect to the tensor and 
scalar fields, respectively, we obtain the full field equa- 
tions of the generic theory described above: 



8ttT* - 



duj{4>) 



^ White-dwarfs typically have sensitivities Sa ~ 10 neutron 
stars have sensitivities Sa ~ 0.2, and black holes have Sa = 1/2; 



see 12711 for detailed calculations. 
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where we have defined T* = T — 20|^ and Dg is the 
curved space d'Alembertian, defined by 

U,^{-g)---d.{{-g)'^a^^d,). (14) 
A detailed derivation of this result can be found in 



C. Massive Brans-Dicke theory: The field 
equations and their wealc-fleld limit 

As we recalled earlier, the effects of a generic coupling 
function on the dynamics of compact binaries have been 
studied fairly extensively by Damour and Esposito-Farese 
Here we are primarily interested in the effects 
of a nonzero mass of the scalar field. In the limit where 
TOs — )■ 0, our final result for the dipolar and quadrupo- 
lar flux can be shown to match Eq. (6.40) in [2l| (the 
monopole contribution vanishes for circular orbits). 

It would be interesting to study a theory with generic 
functional forms for both a;(0) and M((/)), but for sim- 
plicity here we will consider a constant coupling function: 
a;(0) = wbd = constant, as in the usual Brans-Dicke the- 
ory [2^. The scalar field equation then reduces to 



d4> 



2M(0) 
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8ttT* 
3 -I- 2a;BD 



(15) 



In order to get a handle on the effects of the cosmo- 
logical function M(0), let us expand the metric about a 
Minkowski background t]^^ and the scalar field around its 
(cosmologically determined) constant background value 
^0- Following closely the method of [lO], we define small 



perturbations ip, h^^ and 6^1, such that 
(/) = 00 + V?, 

2 \(po/ 
Let us also expand M{(f>) in a Taylor scries about (j)Q: 

Af (0) = Af (,/)o) + M'{<Po)^ + hl"{(t>o)ip^ + ... (17) 

We require that the expanded field equations are con- 
sistent at all orders in (u/c)". Substituting the weak 
field perturbations ((T6)) into the field equations and 
(|14p and examining the leading-order terms under the as- 
sumption of asymptotic fiatness, we find that Af(0o) = 
Af'(0o) = 0. We are therefore left with the quadratic 
term, that endows the scalar field with mass. To see this, 
let us substitute M((j)) = ^M" {(j)Q)(p'^ into the scalar field 
equation, yielding 

- ml{(f) - (j)o) 



8ttT* 



(18) 



(19) 



3 + 2a;BD ' 

where we have defined the mass of the scalar field 

We will see shortly that nis is precisely the parameter 
appearing in Yukawa-type corrections ~ to the 

Newtonian gravitational potential, as well as the ordinary 
mass parameter in the Klein-Gordon equation. Since the 
scalar field is expected to be small, we will neglect cubic 
and higher-order terms in M{(j)), that would introduce 
additional nonlinearities into the scalar field equation. 

In summary, with our choice of coupling and cosmolog- 
ical functions, the field equations of the massive Brans- 
Dicke theory read 



3 -I- 2iUbD 2/ , , \2 
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V,-5^,ng0), (20) 
(21) 



D. The weak- field limit 

Let us use the weak-field perturbations (|T6)) to obtain 
the field equations in the weak-field limit. Expanding the 
left hand side of and imposing the harmonic gauge 
condition 6'^'^ = we find 

R^.-l:9^.R = -1:^,0^.. + ^ - V^.0J^) , (22) 

/ Z 00 ^00^ 



where is the fiat-space d'Alembertian, and we ne- 
glected quadratic and higher-order terms. The tensor 
field equation can hence be written as 

= -leTTT^'^ , (23) 

where r^i/ = 0(^"^T^i/ -I- t^i^. We have collected the 
quadratic and higher-order terms in the perturbations ip 
and 0^1, into the gravitational stress-energy pseudotcnsor 
tfj^^. By virtue of the gauge condition on 6'^'^, we have the 
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useful result that 



r^" = . 



(24) 



Following a similar procedure for the scalar field equa- 
tion, we expand Og(t> in the weak-field perturbations: 



2 00 00 00 00 



(25) 



Substituting this back into the scalar field equation we 
find, as anticipated, the standard Klein-Gordon equation 



where we have defined the source 5* as 



dT\ / 1 



(26) 



5 = - (6 + 4cjbd)"MT-20— 1 - -e- 



f67r* 



/>0 V,a0'" 



n^s0o 



+ 0(6l^6lV,6'<P ,^ ) 



(27) 



II. POST-NEWTONIAN EXPANSION OF THE 
MASSIVE BRANS-DICKE THEORY 

We will now perform a post-Newtonian expansion of 
the scalar and tensor fields. This will allow us to derive 
the Shapiro time delay (section HlH) . the equations of mo- 
tion and periastron shift of compact binaries (section IIVP 
the Nordtvedt effect ( section |V|, and will be required for 
the derivation the period derivative due to gravitational 
radiation (section IVT)) . Before we proceed, it will be con- 
venient to define some auxiliary combinations containing 
wbd that show up repeatedly throughout the calculation: 



1 



2 + wbd ' 
f + ujbb 

2 + ujbd ' 

1 

3 -f- 2WBD 



(28) 
(29) 
(30) 



Furthermore, for our choice of units the cosmologically 
imposed 0o is given by 



4 + 2w 



BD 



3 + 2wbd 



(31) 



Following very closely the method described in [20| (see 
Appendix |X] for details) we obtain 



= ? E — (1 - 2.,)e— +eY: ^^(^a + 2< - 2.^) X (1 - 2s,)e 



2 \ - malTLb 



\e V X (1 - 2sJe-™=- (1 ~ 2.,)e— ^ „ ^ V ^^(1 - 2sje— "00-^ x (l + a(l - 2sb)e— -) 

2 ^ ran TaTab ^ J 

a,o a^o 

^Ea^(^)(l-2.^^) + 0(6)^ (32) 



^^2^— — (l-2s,)e 



ra ' ' 2 

a 

goo = -1 + 200 1 ^ (l + a(l - 2s, 



20-2 ^ (i + a{l - 2sJe-™='-») x (l + a(l - 2s,)e-™='-^) 



, TOgmb r 



1 + ae 



a=tb 

e{sa + 2s',-2sl){l-2sb)e- 



■) X (l + a(l - 2sfc)e-"^''-) - C0o ^^.(e-"-'''" + (1 - 25^)6"™='--) 



E fl + 27 + ^s.e-'"^'-" -t- ^e(l - e-™='-") 



+ E™aC(l-2Sa) X — 



a2 /(2 + 77i,ra)(l-e-™='^")-2m,ra 



2mlra 



0(6), 



X -^-:—ira + 2a(l - 2sa) — 



a " a 

ffy = % + 20o-i*y ^ - a(l - 2sJe-"='-» ) + 0(4) 



m^ra - 1 



0(5), 



(33) 

(34) 
(35) 



a 



I 

In the limit wis 0, the above results reduce to those obtained in the massless Brans-Dicke case [9|. 
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Substituting these results into (|27p we find an cxpres- sion for tlie source S in the near zone, to 0{2): 



2si 



mb{l - 2sb) _j 

n 



III. SHAPIRO TIME DELAY 

Using the post-Newtonian expansion of the metric, we 
can derive an expression for the Shapiro time delay of 
a light ray passing near a massive body. We note first 
that the parametrized post-Newtonian (PPN) formalism 
is not viable when dealing with theories that contain mas- 
sive fields. In fact, Newtonian order terms are modified 
by the presence of massive fields, in the sense that the 
Newtonian potential acquires a Yukawa-like correction of 
the form 



1 



X — x'l 



)d-^x'. (36) 



The impact of this fact for our current purpose is signifi- 
cant: the above potential cannot be expanded in powers 
of 1/r, and the coefficients of modified post-Newtonian 
potentials in the post-Newtonian metric are not con- 
stants, but they have a spatial dependence. Nonetheless, 
we can use the derived metric to obtain an expression for 
the equations of motion of a photon, and use this to ob- 
tain an expression for the Shapiro delay. We will follow 
closely the method described in [2^. A similar calcula- 
tion was carried out by Perivolaropoulos [Tsj : he used a 
different definition of the mass of the scalar field, but his 
results are consistent with those derived here. 
For a photon traveling along a null geodesic. 



0. 



(37) 



To requisite order, 0(2), the equation of motion can be 
written as 



(% + /i|f )uV =0, (38) 



where /i^"^ is the 0{n) order correction to the metric. 
Specializing to a single spherically symmetric source of 
mass M (and negligible sensitivity) at the origin, the 
post-Newtonian corrections to the metric are (from equa- 
tions (031) and ((55)) ) 



Substituting these into 
the photon now reads 



^2U, 



'■)(5., =2i7(l-ae-"»'')(S,,. 

(39) 

, the equation of motion for 
ae-™='')C/) |u|2 =0. (40) 



The unperturbed Newtonian trajectory of the photon will 
simply be x*(i) — xl + n'-{t — tp), where the photon is 
emitted from Xe in direction n at time tg. Let us now 
parametrize the post-Newtonian correction to the tra- 
jectory by Xpjj {t) , where the corrected trajectory is then 
given by x'^{t) = xl + n^{t — te) + Xp^{t). Substitut- 
ing this into the above, we find that the post-Newtonian 
correction to the trajectory satisfies 



dxpp, 



dx 



I 

PN 



-2U. 



dt dt 

Integrating with respect to time, we obtain 



-2 I Udt' . 



(41) 



(42) 



The time taken for the photon to travel from Xe to some 
other point x and back again is hence given by 



At = 2|x- 



Udt' . 



(43) 



The travel time correction 5t due to the Shapiro delay 
corresponds to the second term on the right-hand side. 
Performing the integration, we find for the Shapiro delay 
term 



5t = AM In 



(re + Te • n)(rp - • n) 



(44) 



where the photon is emitted from v^, in direction n, trav- 
els to Tp and back again, M is the mass of the body 
causing the time-delay and rf, is the impact parameter of 
the photon with respect to the source. The mass appear- 
ing in (j44p is not a measurable quantity; what is actually 
measured is the Keplerian mass A/k = A/(l + ae~™=''), 
where r should be thought of as a fixed quantity which 
depends on how the Keplerian mass of the body was de- 
termined. In terms of A/k we have 



8t^ In 
1 -f- ae-'"-'' 

= 2(1 + 7)A/Kln 



(re + Te • n)(rp - • n) 



(re + re • n)(rp - rp • n) 



where in the second line we have defined 



ae 



, (45) 



(46) 
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In the case of the solar system, the r appearing in the 
definition of 7 should be set to lAU, since this is the 
scale associated with the determination of the Keplerian 
mass of the Sun. In any metric theory of gravity where 
the PPN formalism can be applied in a straightforward 
manner, the obtained expression for the Shapiro delay is 
identical to (|^5|) . only with 7 replaced by the PPN param- 
eter 7 (see for example [l^])- We can therefore compare 
7 directly with the observational constraints on 7 from 
Shapiro time delay measurements to obtain an exclusion 
region in the (ojbdi ™s)-plaiie. In section [VII Bl we will 
do precisely this, comparing the derived expression for 
7 to the constraints on 7 from time-delay measurements 
obtained by the Cassini mission. 

Note that in the limit where -^00^ 7 ~^ Ij i-e. the 
GR value of the PPN parameter 7. In the limit where 
TOs — > we have instead 7 — > 7 = (1 + lj-qd) / {2 + wbd), 
i.e. the value of 7 in the massless Brans-Dicke theory. 



and 



=l-^e[2-(l-25a)(l-2s6)e-™='^- 

-(l-2,5,)(l-2s,)e-"='^-] 

+ ic'[l-(l-2s,)(l-2s,)e-™='-- 

-(l-2sa)(l-2s,)e-"=^- 
+ (l-4(s,+,sl-s2)) 

X (l-2sb)(l-2se)e' 



(50) 



Now let us now specialize to a two-body system with the 
center of mass at the origin; to this end let us define 



r = r2 - ri , m = mi + m2 , 
TOim2 



dm = 7712 — mi J M 



m 



(51) 



IV. EQUATIONS OF MOTION AND 
PERIASTRON ADVANCE 



We also write Q12 = G and B12 = B. With this spe- 
cialization made, the equations of motion are found to 
be 



Armed with the post-Newtonian expansion of the 
fields, we are now in a position to obtain the EIH equa- 
tions of motion. From ([7]) , the matter Lagrangian for the 
ath body in the system is given by 



ma{<P) {-900 - 2goiwl - gijvlv'i) ■ 



(47) 



To obtain an n-body action we follow the procedure de- 
tailed after Eq. (11.90) of We substitute the post- 
Newtonian expressions for the metric and scalar fields 
obtained in the previous section and use the expansion 
of ma{4>) in ([TU| . We first make the gravitational terms 
in La manifestly symmetric under interchange of all pairs 
of particles, then wc take one of each such term generated 
in La, and sum over a. To 0(4) we find 



^EIH = 



Gab + SBabvl - ^ Vabc — 



Tab 

a^b 

^{Gab + 6i3ab)Va ' Vfc - ^Gabi'^a ' nab){vb ' Hab) 

(48)" 



g-3gB--l{G-3By 



mr 
— 

1 . -7=, ^ -'^ /-n -rS \ "^2 
77(2^211 + 2^211 ) 7: (2-'l22 + i-'l22 ) 

2 \ ' r 2^ ' r 



2GG'- + {G + 2G) ^v' ^1{AG-G)^{^- nf 



r 

m(r • v)v 



G + 3B+{G-3G)- 
m 



(52) 



where 



G = l- -C[l - (1 - 2si)(l - 252)(1 + m,r)e-™^^n , 

B = i(27 + 1) + - (1 - 2.si)(l - 2s2)(l + m,r)e~ 
6 6 



(53) 



and 

A22 =1 - e[l - (1 - 2si)(l - 2S2)(1 + m,r)e-"^=n 

+ [1 - 2(1 - 2si)(l - 2s2)(l + m,r)e-™='- 

+ (1 - 4(.si + s[ - s?))(l - 2.S2)'(1 + 2m,r)e-2™='-' 

(54) 



where we have defined 



Periastron advance 



nab 



rab 



Gab = l- ^^[1 - (1 - 2.sJ(l - 2sfc)e-'"='--] , 

Bab = ^(27 + 1) + k[l - (1 - 2sJ(l - 2si)e 
3 6 



(49) 



With the equation of motion in hand, we can view 
the post-Newtonian corrections together with the scalar 
Yukawa-like terms as perturbations of the Keplerian orbit 
and employ the method of osculating elements [l^] to 
obtain an expression for the periastron advance of the 
binary system. In contrast to the massless Brans-Dicke 
case (treated in [1^), the integrals that appear in this 
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perturbation expansion cannot be written in closed form, 
so an expansion in powers of the eccentricity e is required 
to obtain closed-form expressions. Fortunately for our 
current purposes we will only require the result in the 
two limiting cases of very light and very massive scalars. 
In the former limit msr <C 1, the periastron advance 
reduces to the masslcss Brans-Dicke result I20ll 



67rm 1 
w = — TT^VQ 



a(l - e2)P 



(55) 



where a and e are the semi-major axis and eccentricity, 
P is the period and V is given by 

v = gB + -g^-- 56 

DO TO 



To Newtonian order, the n-body Lagrangian (l48l) is 
given by 

iEIH = - E TO,(1 - \vl) + I Y: ^Ga,, (58) 

and the n-body equations of motion arc hence 



(59) 



Sb)msrabe ™='''"'fah- 



The relative acceleration of two bodies A and B in the 
field of a third body C is then 



In the limit of a very massive scalar to^t 3> 1 the expres- 
sion for the periastron advance reduces instead to the 
familiar GR result: 



a(l — e/)F 



V. NORDTVEDT EFFECT 

Scalar-tensor theories of gravity predict that massive 
bodies with a significant amount of gravitational self- 
energy do not follow geodesies of the background met- 
ric; in fact, massive bodies with different gravitational 
self-energies will follow different trajectories, leading to 
direct violation of the strong equivalence principle. This 
is known as the Nordtvedt effect, and leads to detectable 
effects in the Solar System. Most notably, it leads to a po- 
larization of the Moon's orbit around the Earth [lO, [l^ , 
which can be constrained using lunar ranging experi- 
ments. Let us look at how this effect arises in the massive 
Brans-Dicke theory. 

The effect is usually parametrized by the Nordtvedt 
parameter ?7n , which can be determined directly from the 
PPN metric of a given theory, and it turns out to be some 
simple combination of PPN parameters. However, as we 
have seen previously, in the case of the massive Brans- 
Dicke theory the PPN formalism is not directly appli- 
cable. However we can extract an "effective" Nordtvedt 
parameter from the equations of motion. To do this, let 
us consider the relative acceleration of a pair of bodies A 
and S, aAB = sa — a-B, in the field of a third body C, 
with rAB TAG and tag — ''bc- The Nordtvedt effect 
will result in an anomalous difference in the accelerations 
of A and B towards C, proportional to the difference in 
the specific grav itational self-energies of the two bodies 
A and B f20ll29l. |30| . Since the sensitivity Sa of a body is 
related to its gravitational self-energy Q.a by Sa = fla/rria 
(in the weak field limit), the extra term arising in aAB 
due to the Nordtvedt effect will be proportional to the 
the difference in sensitivities S = sb — sa- 



SLAB = aB - aA 



(60) 

^ab(toa -f- tob) ^ GBcmc^ , GAcmc. 

2 ""AB 2 ''BG H 2 ^'^C 



' AB 



' BG 



' AG 



C(l - 2sb)(1 - 2sc)m,rBGe-'"-'^^^fBG 



+ i f ^^^^±^e(l - 2sa)(1 - 2sB)TO,e-™='--B I TAB 

_ ^ ™c 
2 ''bg 

+ - 2sa)(1 - 2sc)TO,rAGe-™='-^°fAG ■ 

^ '"AG 

Regrouping terms together appropriately and assuming 
that TAB ^ ''AG, '"AG — ''bGj we can rewrite this as 



aAB 



*rAB I J_ TOGrAG _ TOcTBC \ 



'AB '?^0 V ''AC ''Ig / 

+ [ai - 2.sg)(1 + TO,rAG)e-"='-*°] {sb - sa)^ 



(61) 

rnGfAG 



'AG 



= \ I — I + rj^b — ^ 



' AB 



' AC 



' BG 



'AC 



where the first term is the Newtonian acceleration be- 
tween the two bodies, the second term is the tidal cor- 
rection to the orbit of the system (A, B) and the final 
term (proportional to 5 = sb — sa) is the difference in 
the accelerations of A and B towards the third body C 
due to the Nordtvedt effect (cf. [1^). In the second line 
we have rewritten the third term in the conventional form 
from which the Nordtvedt parameter is usually defined; 
we can then simply read off the effective Nordtvedt pa- 
rameter 



m = C(l + 'msr){l - 2sc)e~ 



(62) 



where r is now taken to be the distance from C to the 
system {A, B). Note that if the Sun were replaced by a 
black hole (sc = 1/2), there would be no Earth-Moon 
Nordtvedt effect. In section I VII CI we will compare the 
effective t^n to the measured value of the Nordtvedt pa- 
rameter provided by Lunar Laser Ranging experiments 
to obtain bounds on (wbd, ^s)- 
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VI. GRAVITATIONAL RADIATION FROM 
COMPACT BINARIES 

A. Tensor radiation 

In this section we will follow very closely the general 
method described in [23| ■ The power radiated in gravita- 
tional waves due to tensor radiation in the Brans-Dicke 
theory is given by 



E 



'327r' 



(63) 



where the angular brackets represent an average over one 
orbital period and Ol^r^ is the transverse-traceless (TT) 
part of 0*^ . 

In order to obtain a formal solution to the linearized 
tensor wave equation (j23p . wc simply fold the source 
T^y with the retarded Green's function of the flat-space 
d'Alembertian operator 



G(i-i',R-r') 



5{t-t' -\R-r' 
|R-r'| 



with the result 



IR-r'l 



(64) 



(65) 



Here the integral over t' has been carried out immedi- 
ately, and the spatial integration region ^ is over the 
near zone. If wc make the assumption that the field point 
in is the radiation zone, such that |r'| ^ |R|, and make 
the slow-motion approximation, we can expand the r' 
dependence of the integrand and write 



^,.^1 - 1 a 



R ^ m! dt^ 

?n— 



■ / T^'^(t-i?,r')(n-r')"d^r', 



(66) 



where n = H/R, and the integration is now over Ai, 
which is the intersection of the world tube of the near 
zone with the constant time hypersurface tj^ ^ t — R 
[3l| . For the purpose of obtaining the power loss due to 
gravitational radiation, we are ultimately interested in 
O'^'^Q. Due to our choice of gauge 9^^q = 6*^0 = 0, and 
hence we only require the spatial components 6^^ , which 
are given (to leading order) by 



4 

= R 



r°°(t-i?,r')/V'^dV. 



2f_ 

Rdt^ 



(67) 



Here we have written the monopole moment of r'-' as the 
time derivative of the quadrupole moment of t^'^ , by ex- 
ploiting the conservation law t''^ together with the slow- 
motion approximation. There can be no contribution 
from the dipolc moment of r"" in ([57]) to order O(^), 



since the time derivative ^ ^ v. The quadrupole mo- 
ment of r'-' only comes in at higher order, and hence we 
only require the leading-order contribution from t'^^: 



Substituting this into ()67p we obtain 



(68) 



(69) 



Specializing to a two-body system with the center of mass 
at the origin using (|5ip . we obtain to the requisite order 



(70) 



where we have used ([5^ to replace P (to leading order) 
where necessary. 

We now need to project ([70)) onto the TT gauge by 
applying the projector 

A(n)ij,fc/ = SikSji - -SijSki - UiUkSji 



1 



1 



1 



+ -jTLkUiSij + -n^rijSki + -UinjUkni , (71) 
which satisfies Aij^ki^kLnm = Aij,nm 111 to O''^: 

0¥t = A(n),,- fci^^' . (72) 

The result is 

R^ 
'327r' 



(73) 



We now note that the only fi dependence in the integrand 
of (j73p in contained in the h-ijM- Performing the integral 
over the solid angle we find 

27r 

K.j^kidn = — {lUikSji ~ 4%4; + ki^jk) , (74) 
where we have used the identity 

n'^ • • • n'- dil = (^21 + ' ^^^^ 

Substituting this result back into ([75)1 we obtain 

E = -|^'^o§(l2 e%9'[o - 4 9\,0\,) . (76) 

At this point we will specialize to a circular orbit, which 
we will parametrize by 



ri = r cos(uj{t — i?)) , 
r2 — r sin(cLi(t — i?)) , 
n = 0, 

Vi — —V sin(a;(i — i?)) , 
V2 = V cos(w(t — i?)) , 
^^3 =0, 



(77) 



(78) 
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where w is the orbital frequency. In addition, let us sup- 
pose that the mass of the scalar is either sufficiently large 
or sufficiently small that variations of G over an orbital 
period can be neglected. Then G will reduce to the mass- 
less Brans-Dicke value in the limit of a low mass scalar Q , 
or to the GR value in the limit of a very massive scalar. 
With these two approximations made, we perform the 
average over one period and obtain the final result for 
the power emitted in tensor gravitational waves in the 
Brans-Dicke theory: 



15 



(79) 



We stress again that this result is only valid in the 
limit where rris is such that either e~™='' ~ 1, in which 
case G reduces to the massless Brans-Dicke value or 
^-m^r _^ jj-^ which case G reduces to the GR value. 



B. Scalar radiation 

The general expression for the radiated power due to 
scalar radiation in Brans-Dicke theory is [20| 



Using the relation {P / P) = —^{E/E) as well as the 
Newtonian result (following from the virial theorem) that 
E = T + V ~ —^nv^ to eliminate v, we finally obtain 
the fractional period decay due to the emission of tensor 
gravitational radiation 



P 
P 



-^^(12-60. 



(80) 



E^^^cl)^\4ojBB+e){(l)V,oV,odn), (81) 



where the angular brackets represent the average over 
one orbital period. 

We can solve Eq. ([26| by using the retarded Green's 
function for the massive wave operator □ — m^: 



G(t-t',R-r') 



6{t - 1' -\R-r'\) 
IR-r'l 



e(t-i'- |R-r'|)- 



v/(t-i')'-|R-r?) 



(82) 



r 



where Ji is the Bessel function of the first kind, and 8 is of this result). Now we can write the general solution to 
the Heaviside function (see [s^ for a detailed derivation (|26|) a,s (p = ipB + Vm, where 



V3i3(i,R) =4 



Sit',r')5{t-t'-\K~v'\,,,^, 



.A' 



|R-r' 



-(Pv'di! , 



.y{t~t'y - |R-r'|2 



.^Mz)S{t-J\R-r'\^ + {^J^,r') 

-4 / (fr' X / ^ : dz . 

Jo 



R-r'P + (-f-)2 



(83) 



Q{t-t' - IR-r'Dd^r'di' 



r 



the spatial integration is over the near zone and 

in the last line we have made the substitution z ~ 
msy/{t-t')^-\R-r'\^. 

Taking the field point to be in the radiation zone 

I 



(|R| 3> |r'|) and making the slow-motion approximation, 
we can expand the r' dependence of the integrand and 
write the general solutions (|83|) as 



4^19" 



R ^ m\ dV 

m=0 



4^19"' 
R ^ m! at" 

m=0 



■J d^r'S{t^ R,r'){n-r' 



d^r'in ■ r')" x / dz 



S{t- ^R^ + {^f,v')Ji{z) 



(84) 
(85) 
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We are now in a position to substitute the post- 
Newtonian expression for the source S into (|M)) and ([55)1 
and obtain an expression for the gravitational waveform 
ip(t, R) in the far-field, slow-motion limit. We must first 
specialize to a two-body system with the center of mass 
at the origin, using (|51|) . Performing the integration and 

I 



retaining terms up to order 0(^^^) and 0(^t) in the 
monopole (m = 0) and quadrupole (m = 2) terms, and 
O(^) in the dipole terms (rn ~ 1), wc obtain (mod- 
ulo time-independent terms that arc uninteresting, as we 
ultimately require ip^ in order to calculate the radiated 
power) 



r(n • v)2 - iri,2 - grm^^!-^ 



(2-Or'--(2A-(2-e)r')-e- 



iPrn = -2ai?- V(r/3 [(n • v)2] - -r/i [v^] - r/3 [^r 



(n • r)^ 



1 

(2A-(2-e)r')/i[- 



25/2 [(n • v)] 



25(n-v) , (86) 



(87) 



Here we have defined 

r' = 1 - si - S2 , 

A EE gr' - e((l - 2si)4 + (1 - 2s2)s'i) , (88) 
and the terms /„ [/(i)] represent the integrals 



where the integration over z has yet to be performed, and 
it is understood that the time-dependent terms in (j86p 
and ([87]) (replacing f{t) in (|89| ) are the components of r 
and V. As in the calculation of the tensor component, we 
assume that Q is approximately constant over an orbital 



period {G — > G), and we specialize to a circular orbit 
parametrized by (j77p . Taking the partial time derivative 
of ip we find 



',0 =2aR~^Gm^l\2Si^—~I2[—] 
^/n n-' 



4r 



(90) 



where we have used ([5^ to replace P (to leading order) 
where necessary. The first and second terms represent the 
dipole and quadrupole contributions respectively; note 
that there is no monopole contribution to leading order 
in the circular orbit case. Substituting this into ()8ip and 
performing the integration over the solid angle via the 
identity ([75]) . wc obtain 



(91) 



-5^(1 - 2Z2(i?; m„ w) + W2{R] m,, w)) + jf^v^^ - '^Z^{R- m„ 2w) + W^iR- m„2uj)) 



I 

where in the second line we have performed the average over one orbital period and we have defined 
I 



Zn{R', nis^uj) = cos{ujR)Cn{R', ms,uj) + sin(cj i?)S'„(i?; ms,ui) 
W„(i?;m,,a;) = |C„(i?; Ws, w)P + |S'„(i?; to,, w)|2 , 



C„(i?; ■ms,Lo) 



r cos(lorJi + {^)A ^dz 

Jo V V m^R' J (I + (^)2)2 



Sn{R',rns^uj) = J sui(^u!R^ 1 + 



^_ Mz) 



-dz . 



(92) 
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To get the total power radiated we must perform the 
integrals in the limit i? — >■ oo in which they have closed 
form solutions. The evaluation of these integrals is dis- 
cussed in Appendix |B] Performing the integrals, we ob- 
tain 



E = - 



15 



4^2 



Q{uj — rris) 



Q{2lO — TTls) 

(93) 



Using again (P/P) 



-liE/E), and E = 



we can eliminate v and find for the fractional 



period derivative due to scalar radiation: 



P 
P 



96 g^fim^ 



12 



Q(2uj — rus) 



'-Q{uj — nis) ■ 



(94) 



Combining this with the result for the tensor gravita- 
tional radiation contribution (j80p . we finally obtain the 
result quoted in Eq. ^ of the introduction. 



VII. OBTAINING BOUNDS ON (cjED m^) 

In this section we compare our results for the period 
derivative of compact binaries, the Shapiro delay and 
the Nordtvedt effect against recent observational data to 
draw exclusion plots in the two-dimensional parameter 
space of the theory, (ms, lobd)- Figure [1] in the introduc- 
tion summarizes our main results. 



A. Bounds from P in compact binaries 

Due to the presence of the difference in sensitivities 
(iS = si — S2) in the dipole contribution to the period 
decay H]), the best candidate systems for drawing ex- 
clusion plots in the (ajBD?™^) plane are mixed binaries. 
White dwarf-neutron star (WD-NS) binaries are partic- 
ularly suitable due to the large difference in sensitivities 
(~ 10-'' and - 0.2 for WDs and NSs, respectively H^). 
To our knowledge, there are three such systems for which 
accurate measurements of P and the other necessary pa- 
rameters have been made (to date): PSRs J0751-I-1807, 
J1012-I-5307, J1141-6545. A summary of the observations 
and the relevant references are provided in Appendix [U] 

In principle, there are two more systems that are of in- 
terest to our current purposes: PSR J1738-I-0333 and 
PSR J1802-2124 [Hj. Both of these systems have very 
small eccentricities, which means that the result derived 
in this paper can be used "out of the box" , with no need 
to generalize our calculation to eccentric binaries. In the 
case of PSR J1738-I-0333, a relatively accurate measure- 
ment of P has been achieved (with error ~ 30%), but the 
masses of the component stars have yet to be determined 



j34{ . PSR J1802-2124 is in precisely the opposite situa- 
tion: the masses of the components have been measured 
to reasonable precision, but a precise measurement of P 
has yet to be achieved. This is anticipated in the near 
future m. 

The general approach to obtaining bounds on 
(uj-BDi'm-s) using observations of the period derivative of 
mixed binaries is as follows. Firstly, we need to write ([T]) 
in terms of the observablcs relevant to the system under 
inspection. For circular binaries, the relevant observablcs 
are the stellar masses (including the mass ratio q) and the 
period. Recasting ([T]) into these observablcs we obtain 



P = P 



GR 



4 Ki 5 _2/27r\ 3 2 

g 3 m 3 — 2S kd 

12 96 \P J 



where Pgr is the prediction from GR, given by 



Pgr = - 



1927r q 1713 

5 {1 + qy 



, (95) 



(96) 



For mildly eccentric binaries, provided the eccentricity 
is small enough, we can get approximate bounds using 
the results obtained here for the circular case. In these 
instances, we can use the measured periastron shift w, 
period, and the mass ratio q of the binary. Recasting ([T]) 
in terms of these observablcs (using the results for the 
periastron advance quoted in section IIV Al and Kepler's 
third law to eliminate m and r) we obtain 



P 



GR 



Ki 



5 2tt .^Kp 
12 ' 16Pw 2 



where 



Pgr 



AqP 



8 



{l + qY 15V3 



(97) 



(98) 



With the predicted P written in terms of the relevant 
set of parameters, we are in a position to compare it 
to observations; the predicted P and observed Pobs are 
consistent to na confidence provided that 



P(e,7n,)| < na. 



(99) 



where a is the combined uncertainty of Pobs and the pre- 
dicted P, and where we should remember that the latter 
is uncertain due to uncertainties in the observablcs (such 
as stellar masses and period). In order to obtain an up- 
per bound on ^ (and hence a lower bound on cjbd) to 
95% confidence for a range of scalar masses, in Figure [1] 
we simply plot the contour in the (ojbd, '7i,s)-plane asso- 
ciated with |Pobs — P(^, "^s)| = 2ct. 



1. Bounds from neutron star-neutron star binaries 

The presence of dipole radiation in mixed binaries sug- 
gests that these should be the best candidates for obtain- 
ing the most stringent bounds on (wbd, ^s)- However, it 
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is worth looking into the bounds that could be obtained 
from observations of neutron star-neutron star (NS-NS) 
binaries. Since the sensitivities of the two component 
stars are nearly identical in this case {S ~ 0: cf. j27j'). 
the expression for the period derivative reduces to 



GR 



12 



(100) 



Expanding to linear order in ^ wc can write this as 



P = P< 



GR 



where 



^ 12 



+ , 



e(2cj-m,)- - 
6 



(101) 



1 



+ -(l-2si)(l-2s2). 



(102) 



We can also write the observed period derivative as 

Fobs ^ PcRil + S) , (103) 

where 5 is the fractional deviation of the observed value 
from the GR prediction. Applying the condition ([M)) 
to the above two equations we find that the predicted 
P is consistent with the observed Pobs to 2a confidence 
provided that 



\^X-S\<2a., 



(104) 



where cr is now the combined uncertainty of x and the 
observed deviation from GR, d. Since the correction x in 
the above is of order unity, we conclude that bounds com- 
petitive with the most stringent bounds obtained here 
(from the Cassini Shapiro delay measurements, and those 
that are expected from a rigorous analysis of PSR J1141- 
6545) could only be obtained from a P measurement to 
a precision of ^ 0.01%. Since the current best measure- 
ments of P for NS-NS binaries are not yet close to this 
precision, we conclude that the bounds that would be ob- 
tained by analyzing such systems would be significantly 
weaker than the most stringent bounds obtained here. 
For NS-NS systems for which P has been measured to 
a precision of ~ 1% (such as PSR J0737-3039 [s^), we 
would expect to obtain relatively weak bounds, compara- 
ble to those obtained here from the quasi-circular WD-NS 
binary PSR J1012+5307. 



B. Bounds from Cassini time-delay data 

The Shapiro delay has been measured in the Solar 
System to remarkable precision by radio tracking of the 
Cassini spacecraft in 2002 [s^]. In theories containing 
only massless fields, these observations are tantamount 
to a measurement of the PPN parameter 7. This has 
been measured to be 



^Cassini 



(2.1±2.3) X 10"'^ = l-t-(5±e. (105) 



As discussed in section IIIIl the mass of the scalar in the 
massive Brans-Dicke theory prohibits us from using the 
PPN formalism in the conventional manner, and the con- 
cept of constant PPN parameters breaks down (see also 
jlSlp. In section IIIIl we derived an expression for the 
Shapiro delay in the massive Brans-Dicke theory, and we 
defined a quantity 7 which is analogous to the PPN pa- 
rameter 7 (at least in the context of Shapiro delay) and 
can be directly compared to the measured value of 7 to 
obtain an exclusion region in the (m^, lobd) plane. Com- 
paring the derived expression for 7 ()46|) with the Cassini 
measurement of 7 (jl05|) . we require that 



2e~5 



a < e" 



{2 + S- 2e) 



(106) 



to 95% confidence. The resulting bounds on ^ and cjrd 
are plotted in Figure [T] by solid black lines (cf. also [l8|). 

Wc find that wbd > 40000 for a range of scalar masses 
nig < 2.5 X 10~^°cV; to 95% confidence. This is around 
one order of magnitude more stringent than the bounds 
provided by the observations of gravitational radiation 
damping in binary systems. In the limit l/m^ <C lAU, 
LdBD can take on any value (as long as wbd > ^3/2). 



C. Bounds from Lunar Laser Ranging observations 

The most precise measurement of the Nordtvedt effect 
to date comes from the Lunar Laser Ranging experiment 

M 

?7^^^ = (0.6±5.2) X 10"^ = ^±e. (107) 

Comparing this observed value to Eq. and neglecting 
the small sensitivity of the Sun, we require that 



+ msr)e~"'^'' - S\ < 2e 



(108) 



to 95% confidence, from which we obtain exclusion re- 
gions in the (mj, wbd) plane. These are displayed in 
Figure [1] by dotted blue lines. 



VIII. CONCLUSIONS 

In this paper we set constraints on massive Brans- 
Dicke (or Bergmann- Wagoner) theories with an action 
of the form ([6]) with Ltj{(t>) = wbd, assuming that only a 
mass term is present in the expansion of M{(p) around 
some cosmologically imposed value (f)Q. In particular we 
computed the orbital period derivative for quasicircular 
binaries. From an observational standpoint it will be 
important to generalize our work to eccentric (and pos- 
sibly spinning) binaries, that could yield more stringent 
constraints on scalar-tensor theories. It will also be in- 
teresting to explore possible bounds on massive scalar- 
tensor theories that could result from Earth- and space- 
based gravitational-wave observations of compact bina- 
ries, along the line of Refs. [39l - l4l| . 
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A second obvious generalization of our work will con- 
sist in relaxing our assumptions on the form of w(0) and 
M((j)). More generic assumptions on these functions are 
necessary for a deeper understanding of binary dynamics 
in the context of modified gravity models that try to ex- 
plain cosmological observations. It will be interesting to 
verify whether time- varying boundary conditions on the 
scalar field may lead to interesting binary dynamics [l6j . 

Last but not least, full numerical relativity simulations 
of compact binaries in scalar-tensor theories are under in- 
vestigation by several groups (see e.g. d^H^). Numer- 
ical progress in evolving binary dynamics in alternative 
theories is important, as it could reveal strong-field effects 
that may be inaccessible to post-Newtonian or perturba- 
tive calculations. 



1. Step 1: ip to order 0(2) 

To the lowest order, the scalar field equation (PT|) re- 
duces to 



(A2) 



Expanding the modified stress-energy tensor T* to lowest 
order wc obtain 



a 



(A3) 



Substituting this into (|A2[) . we find the solution for to 
0(2) 
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2. Step 2: ftoo to 0(2) 

The 00-component of the tensor field equation (|20p to 
0(2) is given by 



^00 = -^V^/ioo = Stt ( Too + ^ ) - ^V^ ( ^ ) . (A5) 



2 J 2 V^o 



In a similar fashion to (jA3p , we write for the stress-energy 
tensor (to lowest order) 

T = -Too = - ^maS^{x - Xa) 

a 



Appendix A: Post-Newtonian expansion of the 
scalar field and of the metric 



Here we provide details of the derivation of the post- 
Newtonian expansions ([5^ . ([M]). ([M)) and (1551) . Wc fol- 
low very closely the method outlined in [20|. For our 
current purposes, we must solve the field equations (|20)) 
and (PT|) to the following orders: 



^~0(2) + 0(4), 
hoo-0{2) + 0{4.), 
h^r ^ 0(3), 
h,j ~ 0(2). 



Substituting this along with the derived 0(2) expression 
for Lp into (jASp . we obtain the 0(2) solution for H^q: 



hoo = 2<^-i ^ !^ [1 + a (1 - 2s,) e-™^'^"] . (A7) 



by 



3. Step 3: hi^ to 0(2) 



The ij-component of the Ricci tensor to 0(2) is given 



R^, = - hoo,^, + /it, - h^,kj - hlu^)■ (A8) 



(Al) 

Imposing the gauge condition 



We will do this in a number of steps, as described in the 
following. 



(A9) 



16 



we can write the zj-componcnt of tensor field equation 
pl| to 0(2) as 



(AlO) 



Using Eq. (|A6p for the stress-energy tensor and substitut- 
ing the derived 0(2) expression for tp into the expression 
above, we obtain the solution for hij to 0(2): 

h,, ^ <5,,20o 1 ^ [1 _ a(l - 2sa)e-'"='-»] . (All) 



4. Step 4: hoi to 0(3) 



by 



The Oi-component of the Ricci tensor to 0(3) is given 



Imposing the further gauge eondition 



"2^00.0 



this reduces to 



1 



,0i 



h' 

12 



fe.Oi- 



(A12) 
(A13) 
(A14) 



We can hence write the Oi-component of the tensor field 
equation ([20|) to 0(3) as 



o'To^ + U^] +i:h%n.. (A15) 



'0 / ,0i 



12 



The Oi-component of the stress-energy tensor to lowest 
order is given by 



a 



(A16) 



In order to write the ip^oo term in the form V^x, we must 



find a particular solution to V x 



''/r. Taking 



care to ensure that the chosen solution x is such that the 
correct limit is obtained as rrig 0, we write 



TUsra - 1 



■mlra 



(A17) 



Noting also that V^(ra/2) ~ Tq, we can re- write the sec- 
ond and third terms in (jAlSP as 



5. Step 5: tp to 0(4) 

Expanding to 0(4) and recalling the definition of 
6""" in Eq. (HH), we obtain 



(A19) 



The scalar field equation ((2T|) to 0(4) can hence be writ- 
ten as 



(V^-^^)^^8.a0-r*(l-l.-£ 



Vf ) . (A20) 
Po J ,00 V -Po . 



Expanding the modified stress-energy tensor to the re- 
quired order we find 



T*=^m„[(2s„-l) + -(l-2sa)«2 

a 

-^(l-2,s,)(?- ("2.1-2.^ + 1)^1. (A21) 



In a similar fashion to step 4, we wish to write the (p^oo 
term in the form (V^ ^f^s)Xi require a particular 

solution to (V^ — TOs)x = e~™='"/r such that the correct 
limit is obtained as mg — )■ 0. To this end we write 



(V^-m^) 



1 -e- 



2m, 



(A22) 



and hence write the second term on the right hand side 
of Eq. ((X20I) as 



f_ 

00 / ,00 



= (V^ - m2)e^,7i,(l - 2s,)5oao 



1-e- 



2ms 
(A23) 



^-a(l-2s,)- 



Tflsra - 1 



7- =eV"ia(l-2s„)9,9o 

(A18) 

Substituting these into (|A15[) . we obtain the solution for 
/loi to 0(3) given in Eq. 



The third term on the right hand side of Eq. (|A20p can 
be re-written (to the required order) as 



Substituting the above along with the derived 0(2) ex- 
pressions for iy9, ft-oo s-nd hij into (jA20l) . we obtain the 
result given in Eq. ([5^ . 
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6. Step 6: hoo to 0(4) 

The 00-component of the tensor field equation (|20p to 
0(4) is given by 



1 , 



0/ ,00 



1 w 
-V/iooV^ 

2 <Po 




(A25) 



where we have used the gauge conditions (jA9| and ()A13P 
to reduce the expression for i?oo into a convenient form. 
The term involving the stress-energy tensor on the right 
hand side is given (to the required order) by 



00 



:9ooT 



'2^ 



ma(5^(x- Xa) 



3 2 
— V 

4 " 



4 

(A26) 



Using this along with 
1 



V «go 



00 2 



<y5 



yo 



/iqoV^/iooj 



(A27) 



Eq. (jA25P can be re-written as 
1 

, To 



20^1^ ma 



1 



3 2 
— V 

2 " 



In ^ 
2 (?io 



00 2 00 



00 



2 ^ '^00 



2;iooV2/ioo - 2^W^hoo - I 



00 00 



00 



00 7,00 

(A28) 



Using Eq. (|A17[) to re-write the term involving (p^Qo in 
the form V^x, and substituting in the derived 0(2) ex- 
pressions for hoo, hij and cp, and the 0(4) expression for 
(f, we obtain the result presented in Eq. ([33]) . Note that 
there are two contributions to the term involving the sec- 
ond time derivative in Eq. (|33p : one contribution from ip 
(to 0(4)) and one from (f^o- 



Appendix B: Evaluation of integrals Cn{R\rns,ijj) and 
Sn{R\ma,u)) arising in the derivation of the period 
derivative due to scalar radiation 

In reaching the final expression for the power emitted 
in scalar gravitational radiation (|94p . we were required 



to evaluate the integrals C„(i?;TOs,w) and S'„(i?; TOs, w) 
defined in In this appendix we give details of the 

evaluation of these integrals. 

Since we are interested in the gravitational radiation in 
the far zone, we only need to determine the asymptotic 
behavior of these integrals for i? — > oo. Substituting 
u = ■\/T+(z/m7R)^ into we obtain 



CniR]ms,uj) = rUsR 
Sn{R; ms,Lo) = rUsR 



du 



du 



cos{ll!Ru) JilmsRy/u^ — 1) 
8m{ujRu) Ji[miiR\/v? — 1) 



,n — 1 



1 



(Bl) 



We will discuss the evaluation of C„ only, as the evalu- 
ation of Sn proceeds in exactly the same way. To begin 
with, let us choose some e such that rrisRe 3> 1 while 
LoRe^ <C 1 and split up the integral into an integration 
from 1 to 1 -f e^/2 and from 1 -|- e^/2 to oo. In the first 
integral, as the argument of the cosine is nearly constant 
we can approximate 

^ cosiujRu) Ji{msRVu^ - 1) 

171 a R, I du - 



,n— 1 



1 



« cos(wi?)(l — Jo{msRe)), 
(B2) 

with the zeroth order Bessel function Jo given by its 
asymptotic value 



JoinisRe) 



2 cos(m,,i?e - 7r/4) 
TT ^/nisRe 



(B3) 



For the second integral, we can approximate the Bessel 
function Ji by its asymptotic value 



2 cos(a; - 37r/4) 



Jiix) 



and hence the integral can be approximated by 



— y/nisR i du 



cos{ojRu) cos{msR\/v?' 



(B4) 



37r/4) 



(B5) 



Performing an integration by parts exactly cancels the 
corresponding boundary term in (jB2|) : this is not sur- 
prising, since we expect that the result should not depend 
on the value of e. In analyzing the above integral, then, 
we can neglect all terms arising from the lower endpoint 
(since a full analysis will show that they will exactly can- 
cel the terms arising from the upper endpoint in (jB2[) ). 
Wc are interested in the leading asymptotic behavior of 
the above integral; we hence require the asymptotic be- 
havior of integrals of the type 



du 



+£^2 u«-l^u2 _ I 



(B6) 
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where 



p{u) — iR{niUJU + n2ms\/v? — 1) — m237r/4, (B7) 

with ni.2 ~ ±1- The part of the integration contour 
which gives the dominant contribution is determined by 
p(u) and the relative sizes of a; and m^. Let us deal with 
the two cases a; > to, and cj < to, in turn. 



the branch of \Jv?- — \ that we are using. This will then 
determine the constant phase contours and the location 
of the saddle points. The asymptotic behavior obtained 
in the end will of course be independent of the choice of 
branch. Writing u = w -\- iv and \J v? — \ = J7 + iV ^ we 
define the principal branch of the square root to be 



U + iV w,v>0 

U-iV w>0,v<0 

-U -iV w,v <0 

'U + iV w<0,v>0 



(BIO) 



For Til 



-712, p{u) has a stationary point at a 



Ltj I \J up- — ml and we can apply the method of stationary 
phase (see e.g. Since only a small region around the 

stationary point contributes to the integral, expanding 
the exponent around a gives the leading-order behavior 



oP(a) 



1 2 I 

— \—\/msR- 

4V^ a^-i^/^iT^" J-S 



+s 



1 / ^uP- — m 

2 V ~ 



2^ n-1 



(B8) 



For n\ = ri2, p(u) no longer has any stationary points 
in the integration domain, so the leading order behavior 
is obtained by integration by parts. Since the integrand 
goes to zero at the upper endpoint -|-cx), the only con- 
tribution will come from the lower endpoint, which as 
we have discussed must exactly cancel with the corre- 
sponding terms from (|B2|) . The complete leading order 
behavior of C„(-R; m,,, a;) (and similarly S„{R; imis, ui)) for 
Lo > TOs is hence given by 



The second branch of the square root is then simply the 
negative of (|B10|) . 

For ni = n2, there are no saddle points in the prin- 
cipal Riemann sheet, so we proceed with a straightfor- 
ward integration by parts. As before, since the integrand 
vanishes at the upper endpoint -|-c» the only contribu- 
tion comes from the lower endpoint, and this will exactly 
cancel the corresponding contribution from (|B2[) . 

The case ni = — n2 is somewhat more challenging. In 
the principal Riemann sheet, we find two saddle points 
on the imaginary axis at u = b± — ±iuj/ y^rn^ — up . At 
these saddle points the imaginary axis intersects another 
constant ph ase conto ur that closes on the real axis at 

U = ±TO, 



— -j cos{R^/i 



S'„(i?; TOs, uj) ^ sin(wi?) — 



2s n-l 



\ UJ 



sin(i?^ — 
(B9) 



saddle point b 



UJ < rUs 



The case uj < nis is somewhat more subtle. Now the 
first derivative of p{u) can only vanish on the imaginary 
axis. We must therefore consider the analytic properties 
of p{u) and use the method of steepest descent [4J| . The 
central idea behind this method is to deform the integra- 
tion contour in such a way that it follows lines of constant 
phase (lines of steepest descent), in the hope that along 
the new contour the integral may be evaluated asymp- 
totically. Firstly, we must take care of the fact that p{u) 
is not analytic in the complex plane, owing to the fact 
that the square root term makes it double- valued. Work- 
ing on a two-sheeted Riemann surface we can still apply 
the method of steepest descent, provided the deformed 
contour does not include either of the branch points that 
appear on the real axis at u = ±1. Since we are extend- 
ing our exponent into the complex plane, we must first fix 



■; +1 



XT" 




FIG. 2. Deformed integration contour C\ + C2 + C->, + Ci along 
lines of steepest descent for p{u) = iR(uju — ma\/v? — 1) -|- 
i37r/4 (corresponding to n\ = -f 1, 712 = —1) on the two sheets 
of the Riemann surface associated with \/u^ — \. Here only 
the saddle point on the negative imaginary axis contributes. 



Let us consider the case n\ 



-1. The original inte- 
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gration contour runs along the real axis from 1 + 6^/2 to 
+00. We now deform the contour by going along constant 
phase lines in the direction in which the real part of the 
exponent p{u) decreases (see figure [5]). Starting from the 
lower endpoint 1 + we follow the contour Ci with 
phase Lu — nise into the lower half of the complex plane, 
and then through the branch cut onto the second sheet 
of the Ricmann surface. On this sheet there arc no sad- 
dle points, and our constant phase contour approaches 
(w — mge)/{mg + w) ± too. We can then connect it to 
the imaginary axis by a path C2 parametrized by i + iT, 
with t running from (w — mse)/(ms + cj) to and con- 
stant T ^ 1. The contribution from this path vanishes as 
T — > 00. The integration contour then follows C3 along 
the imaginary axis through the branch cut onto the first 
sheet, through the saddle point at — —icu / ^ — 
and towards —00. From there it can be closed onto the 
positive real axis by the path C4 (analogous to C2) that 
ultimately makes no contribution. Since the integrand 
vanishes as u — > +00 the upper endpoint of the integra- 
tion is once again unimportant. For the case n\ = — 1, 
we can proceed in a similar fashion, only this time the 
deformed contour will pass through the saddle point 6+. 

We now have all of the ingredients we need to eval- 
uate the integral. Since our deformed contour together 



with the original contour does not include any of the 
branch points, we can still apply the Cauchy theo- 
rem and approximate the integration along the contour 
C*! -I- C2 -I- C3 -I- C4 to obtain the asymptotic behavior of 
/ as i? — 00. The two crucial regions of the deformed 
contour are the lower endpoint on the contour C\ and the 
saddle point on the imaginary axis (&_ for n\ ~ -|-1 or 
6+ for n\ = — 1). Even though the contribution from the 
saddle point is sub-dominant, it does give the asymptotic 
behavior of the original integrals ([M)) that we require; re- 
call again that the contribution from the lower endpoint 
will be exactly canceled by the contribution from (|B2p . 
Integrating through the saddle point h± along the imag- 
inary axis and parameterizing u = i{b± + 1), we obtain 




where I± corresponds to rii = ±1. The complete 
leading-order behavior of C„(i?;TOs,w) (and similarly 
Sn{R] ms,u])) for a; < is then given by 



C„(i?;m.,a.) ^ cos(c.i?) - ( V^^^) 



,2x Tl-1 



Sn{R;m„Lu) - sm{ioR) - (X^^l-^) 



(B12) 



Appendix C: Observational data on compact 
binaries used in this paper 



TABLE I. Parameters relevant to the binary system PSR 
J1012-F5307 [H. 



Period, P (days) 
Period derivative (observed), pobs 
Period derivative (intrinsic), P'"'"^ 
Mass ratio, q 
NS Mass, mi (Mq) 
WD Mass, m2 (Mq) 
Eccentricity, e (10"*^) 



0.60467271355(3) 
5.0(1.4) 10~" 
-1.5(1.5) 10"^-* 
10.5(5) 
1.64(22) 
0.16(2) 
1-2(3) 



The relevant parameters for this system are listed in Ta- 
ble m The parameter values are taken directly from [4^ . 
The mass ratio and individual masses were determined 
in [47| . and the intrinsic period derivative, corrected for 
Doppler effects, was determined in j45| . 



Using the parameters listed in TableUand Eq. ((96)) . the 
value of the period derivative predicted by GR is given 
by 



Pat 



1927r 



q 



5 (l + g)2 



- ■ = -1.1(2) X 10- 



(Cl) 



c. PSR J1012+5307 

PSR J1012-I-5307 is a 5.3ms pulsar in a 14.5hr quasicir- 
cular binary system with a low-mass WD companion j46j . 



Using the method described in section IVII Al we obtain 
the bound on ^ (and hence wbd) as a function of the 
scalar mass which is displayed in Figured] by a solid green 
line. In particular, we find a lower bound wbd > 1250 for 
TOs < 10~^°eV. The limiting factor here is our ability to 
obtain a precise value for the intrinsic period derivative, 
once Doppler effects have been accounted for. 
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d. PSR J0751+1807 



PSR J0751+1807 is a millisecond pulsar in a 6hr cir- 
cular binary system with a helium WD companion psj . 
The period derivative has been measured to ~ 15% pre- 
cision, after kinematic corrections have been made. How- 
ever, the determination of the masses of the stars in this 
system has proved to be more of an issue. Assuming GR 
to be true, Nice et al. [i^ used combined observations 
of the Shapiro delay and orbital period derivative to con- 
strain the masses of the component stars to a precision 
of ^ 10%. Unfortunately, in the context of using the 
measured period derivative to constrain modified theo- 
ries of gravity, we cannot assume GR in the calculation 
of the masses. The solution to this issue is to use only 
the observations of the Shapiro delay to constrain the 
masses, and use these masses in conjunction with the ob- 
served P to compare theory with predictions. The prob- 
lem with this is that using the Shapiro delay alone pro- 
vides a very weak constraint on the masses, with ~ 100% 
uncertainty for each of the two components. Nonethe- 
less, we could perform an analysis similar to that done 
for PSR J1012-I-5307. Given the large uncertainties as- 
sociated with this system, however, we expect that the 
bounds obtained from such an analysis would be very 
weak and would not provide us with any further insight, 
and for this reason we have neglected this system. 



TABLE II. Parameters relevant to the binary system PSR 
J1141-6545 [lOl. 

Period, P (days) 0.1976509593(1) 

Period derivative (observed), P°''" -4.03(25) 10"^^ 

Period derivative (intrinsic), P'""^ -4.01(25) 10"^^ 

Mass ratio, q 1.245(14) 

NS Mass, mi (Mq) 1.27(1) 

WD Mass, m2 (Mq) 1.02(1) 

Eccentricity, e 0.171884(2) 

Periastron advance, Co {°yr^^) 5.3096(4) 



e. PSR JII4I-6545 

PSR J1141-6545 is a 394ms pulsar in a moderately 
eccentric binary system with a WD companion |5l| . The 
relevant parameters for this system are displayed in Table 
im and are taken directly from ISOl. The masses of the 
WD and NS were determined by |50| . 

This system is comfortably the most useful in the con- 
text of putting bounds on (wbd, f^s), and in constraining 
alternative theories of gravity using observations of the 
orbital period derivative in general. P has been mea- 
sured to remarkable precision, currently ^6%, and this 
is expected to improve further to ~2% by 2012 (HOj. The 
other necessary parameters for our purposes, the masses 
and the periastron shift, have also been measured to ex- 
cellent precision, so the total uncertainty in the system is 
(relatively) very small. Unfortunately this system does 
not have negligible eccentricity, so the result for P de- 
rived here does not strictly hold. In order to do a full 
and accurate analysis of this system, the result ^ must 
be generalized to cover eccentric binaries. For the mo- 
ment we present a rather crude analysis of this system 
where we neglect the eccentricity, to find at least a ball 
park estimate of the bounds that we may expect to ob- 
tain once a full analysis is performed. Using the above 
parameters and equation (|M)) . the value of the period 
derivative predicted by GR is given by 



P, 



GR 



Aq 



8 



.( 



P 



[l + qf 1573^2^ 



= -3.440(3) 10"" . 

(C2) 



Using the method described in section fVlI Al we obtain 
the bound on ^ (and hence wbd) displayed in Figure [T] 
by a solid blue line. Once we account for eccentricity 
in a proper way, this system is very likely to provide the 
most stringent bounds among all of the binaries observed 
so far. 
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